%% Description

% This script simulates the transient response of a capacitive voltage
% transformer.

%% Initialization

% Define CVT parameters, including: load resistance R0 [Ohm], resonance
% suppression inductance Lf [H], resistance Rf [Ohm], inductance Lf [H] and
% capacitance Cf [F] of the resonance suppression circuit, resistance R
% [Ohm] and inductance L [H] of the compensation reactor and sum of dividing
% capacitances C [H] 

R0 = 1.03997e5 ;    
Lf = 315.3 ;      
Cf = 0.0285e-6 ;  
Rf = 77379 ;        
R = 3289 ;          
C = 9.1605e-8 ;    
L = 76.136 ;        

% Define the CVT transfer function coefficients.

A3 = Lf*Cf*Rf*R0*C ;
A2 = Lf*R0*C ;
A1 = Rf*R0*C ;
B4 = Lf*Cf*(Rf+R0)*L*C ;
B3 = L*C*Lf + R*C*Lf*Cf*(Rf+R0) + Lf*Cf*Rf*R0*C ;
B2 = L*C*(Rf+R0) + R*C*Lf + Lf*Cf*(Rf+R0) + Lf*R0*C ;
B1 = R*C*(Rf+R0) + Lf + Rf*R0*C ;
B0 = Rf + R0 ;

% Define the CVT transfer function.

num = [A3 A2 A1 1] ;
den = [B4 B3 B2 B1 B0] ;
sys = tf(num,den) ;

% Define the simulation parameters.

f = 60 ;            
w = 2*pi*f ;
T = 1/f ;           
t_fault = 0.05 ;    
t_start = t_fault - 5*T ;    
t = t_start : 1e-5 : 10*T ;

% Define the 3PT1/6PT1 class limits at different times after the fault.

Ts_limits = t_fault + [20e-3 40e-3 60e-3 90e-3] - 0.02 ;
ratio_limits = [10 10 10 10]/100 ;

%% Fault at zero crossing

% Define the primary voltage.

Vp = sin(w*t) ;
Vp(t >= t_fault) = 0 ;

% Simulate for the secondary voltage.

Vs = lsim(sys,Vp,t-t_start) ;

% Plot the results.

figure ;
plot(t,Vs,'LineWidth',2);
hold on ;
plot(t,Vp,'--','LineWidth',2) ; 
plot(Ts_limits([2 end]), [ratio_limits(1), ratio_limits(end)], 'r-.', 'LineWidth', 1.5) ;
plot(Ts_limits([2 end]), -[ratio_limits(1), ratio_limits(end)], 'r-.', 'LineWidth', 1.5) ; 
xlabel('Vrijeme [s]','FontWeight','bold') ;
ylabel('Napon [r.j.]','FontWeight','bold') ;
legend('Primarni napon','Sekundarni napon','Granice tranzijentnog odziva klase T1','Location','northeast') ; 
xlim([1.5*T t_fault+70e-3]) ;
ylim([-1 1]) ; 
grid on ;
set(gca,'FontName','LM Roman 10','FontSize',12) ;
set(gcf,'color','w') ; 

%% Fault at minimum voltage peak

% Define the primary voltage.

phi = -pi/2 - w*t_fault ;
Vp = sin(w*t + phi) ;
Vp(t >= t_fault) = 0 ;

% Simulate for the secondary voltage.

Vs = lsim(sys,Vp,t-t_start) ;

% Plot the results.

figure ;
plot(t,Vs,'LineWidth',2) ;
hold on ;
plot(t,Vp,'--','LineWidth',2) ; 
plot(Ts_limits([2 end]), [ratio_limits(1), ratio_limits(end)], 'r-.', 'LineWidth', 1.5) ;
plot(Ts_limits([2 end]), -[ratio_limits(1), ratio_limits(end)], 'r-.', 'LineWidth', 1.5) ; 
legend('Primarni napon','Sekundarni napon','Granice tranzijentnog odziva klase T1','Location','northeast') ; 
xlabel('Vrijeme [s]','FontWeight','bold') ;
ylabel('Napon [r.j.]','FontWeight','bold') ;
xlim([1.5*T t_fault+70e-3]) ;
ylim([-1 1]) ; 
grid on ;
set(gca,'FontName','LM Roman 10','FontSize',12) ;
set(gcf,'color','w') ; 